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We revisit arguments that the structure function F2 P (x,Q ) for deep inelastic y*p scattering is 
hadronic in nature, similar to vector dominance for real yp scattering. Thus, like all other known 
hadronic scattering, including jp, the growth of F^ (x,Q 2 ) is limited by the Froissart bound at 
high hadronic energies, giving a ln 2 (l/a;) bound as Bjorken x — > 0. The same bound holds for the 
■ individual quark distributions that contribute to F^ v . In earlier work, we presented a very accurate 

CNl ' global fit to the combined HERA data based on a fit function which respects the Froissart bound at 

small x, and is equivalent in its x dependence to the function used successfully to describe all high 
energy hadronic data. We discuss the extrapolation of the fit to the values of x down to x = 1CP 14 
encountered in the calculation of neutrino cross sections at energies up to E v — 10 17 GeV, an 
extrapolation of a factor of ~3 beyond the HERA region in the natural variable ln(l/a;). We show 
O |. how the results can be used to derive the relevant quark distributions. These distributions do not 

D ■ satisfy the "wee parton" condition, that they all converge toward a common distribution xq(x,Q 2 ) 

tIh ' at small x and large Q 2 which is just a multiple of F^, but still give results for the dominant 

neutrino structure function F 2 V ^ which differ only slightly from those obtained assuming that the 
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wee parton limit holds. 

PACS numbers: 12.38.Bx,12.38.-t,13.60.Hb 



INTRODUCTION 



The experimental program at HERA, the electron-proton collider at DESY, probed deep inelastic scattering (DIS) at 
small values of the Bjorken variable x, given in terms of the proton momentum p and the electron momentum transfer 
q in the scattering by x ~ Q 2 /2p ■ q. The measurements covered the range 10 -6 to 10 _1 in x, with a corresponding 
range from 1 GeV 2 to 5000 GeV 2 for the virtuality Q 2 = —q 2 . The results of the extensive measurements by the HI 
H ■ and ZEUS detector groups show that the structure function F^ix, Q 2 ) rises rapidly as x decreases with Q 2 fixed. 

A series of papers over the past few years [H-Q made the case that the reduced cross section in DIS, basically 
the structure function F2 P (x,Q 2 ), is hadronic in nature and satisfies a saturated Froissart bound, implying that 
F2 P (x,Q 2 ) — > constant x ln(l/x) 2 for x — > with Q 2 fixed. The basic argument is that the structure function 
is determined by the total cross section of an off-shell gauge boson 7* on the proton at a y*p center-of-mass energy 
squared W 2 = s, where s = (p + q) 2 is the usual Mandelstam variable, and is thus subject through analyticity and 
unitarity constraints to the saturated Froissart bound on total hadronic cross sections cr(s) — > <7q log 2 s as s — > 00. The 
picture is compelling in light of its success in describing hadron-hadron and photon-hadron total cross sections over 
many orders of magnitude with the same basic functional form [Bj]. Moreover, its predictions for the proton-proton 
and proton-air cross sections at the LHC and the Pierre Auger Observatory Q, respectively, are confirmed by 
these new high energy experiments (]~ol-|l~2| . 
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In this paper, we investigate the implications of this bounded behavior for the ultra-small x, large Q 2 limit of the 
quark distributions in the proton using our recent Froissart-bounded fit to the combined HERA data. We show that 
the individual quark distributions can be derived to reasonable accuracy directly from F^ p ', and present the results 
obtained using the extrapolation of our fit to ultra-small x. The extrapolation should be reliable: the fit function 
becomes a simple quadratic in the natural variable v = ln(l/a;) with well-determined coefficients for x small or v 
large, and the extrapolation necessary to reach x = 1CP 14 (v = 32.2) involves only a factor ^3 increase in v from the 
upper values attained in the HERA region. 

In earlier calculations [H, Q of ultra high energy (UHE) neutrino-nucleon cross sections, it was assumed that the 
quark distributions could be treated in a "wee parton" limit in which the individual quark distributions all converge 
to a common quark distribution xq(x, Q 2 ) at large Q 2 and small x. This allowed the replacement of individual quark 
distributions in the neutrino cross sections by xq, given in terms of by xq = Fj^/^^ef, where the are the 
quark charges and the sum runs over the active quarks and antiquarks. We show here that the wee parton condition is 
not satisfied by the quark distributions determined directly from our fit to F^ p . However, the relations between 

and the corresponding charged- and neutral-current structure functions F%^ and FO^ in neutrino and antineutrino 
scattering which were derived using the wee parton assumption continue to hold to high accuracy. 

In Part II of this work, the companion paper [l3| to this one, we extend our earlier calculations of UHE charged- and 
neutral current neutrino-nucleon cross sections (Jq C (E v ) and a u NC {E u ) up to E v = 10 17 GeV, currently the highest 
energies where there are experimental bounds on cosmic neutrino fluxes [Til , [l5| . These calculations require input at 
large Q 2 (Q 2 > 10 4 ) and small x, down to x ~ 10 -14 . The results presented here are highly significant for that work. 
The results in Part II extend our earlier work both by pushing the neutrino energy higher, to E v = 10 17 GeV, and 
by including the contributions of the b quark. 

The present paper is organized as follows. In Sec. |Hl we develop our arguments with respect to the relevance of 
the Froissart bound and its consequences for the form used in our parameterization of the small- x HERA data for 
F2 P {x, Q 2 )- We then summarize the results of our fit [l6| to the combined HERA data (T3), and list the parameters, 
their errors and the significance of the fit. In Sec. IIII Al we discuss the derivation of the individual small x quark 
distributions from our analytic fit to F^ix, Q 2 ). This requires information on the singlet quark distribution F s {x, Q 2 ) 
which can be expressed in terms of a "bare" structure function F^ p and a set of non-singlet quark distributions. We 
show in an Appendix how F^ p can be related readily to our fit to K7 P , and develop simple results for the effect of 
QCD evolution on the non-singlet distributions in Sec. IIII Bl Our results for quark distributions at ultra-small x are 
given in Sec. IIV Al We discuss their implications with respect to the wee parton picture and neutrino cross sections 
in Sec. IIVB1 where we show that F^^, the dominant neutrino structure function, can be expressed directly in terms 
of FJ P to good approximation despite the failure of the wee parton limit used in earlier discussions of this connection 
[1,3. In Sees. ITV~Cl and HVDl our results are compared with other extensions of quark distributions to ultra small x 
based on conventional parton-level fits to the HERA data. We summarize and draw conclusions in Sec. [Vj 



II. EXTRAPOLATION OF F^ p (x,Q 2 ) TO ULTRA SMALL x 

As we emphasized in the Introduction, for the energies E v of interest for UHE neutrino cross sections, we must know 
F2 P (x, Q 2 ) at values of x many orders of magnitude below the range where it has been measured at HERA. To cover 
the highest energy range reached by neutrino telescope searches [HHHl, F v ~ 10 17 GeV, requires an extrapolation of 
eight orders of magnitude in x below the lowest values x ~ 10 -5 -10 -6 encountered at HERA. While this is actually 
only a factor of ~ 3 increase in the maximum value of the natural variable v = ln(l/x), it is still essential that the 
form used to extrapolate F^ p (x, Q 2 ) be consistent with both the asymptotic limiting behavior expected theoretically 
and the present data. If a fit to the data indicates that the measured structure function is already consistent with 
the limiting asymptotic form, the extrapolation may be expected to be robust; our approach, which we summarize 
here, has this feature. It is the best one can do. 

Up to well understood corrections, the reduced cross section reported experimentally [l?} is equal to F2 P (x,Q 2 ). 
This structure function contains all of the strong interaction dynamics, a point made clearly in Ref. From the 
point of view of vector dominance, expected to hold at high energies and large Q 2 (l8| , it is just the extension of real 
7P scattering with photon 4-momentum squared q 2 = 0, to virtual j*p scattering with a virtual photon 4-momentum 
squared q 2 = —Q 2 < 0. We denote the Mandelstam variables for j*p scattering by s, t, and u. 

There appears to be no obstacle to the continuation from q 2 = to an off-shell q 2 = — Q 2 < 0. The total cross section 
for virtual photon-nucleon (j*p) scattering is proportional to the virtual forward Compton scattering amplitude for 
zero momentum transfer to the nuclcon, t = (p — p') 2 = 0. A complete all-orders analysis of the latter in perturbation 
theory [r| shows that it is real analytic in the Mandelstam variable s = (p + q) 2 for the scattering of a virtual photon 
from a nuclcon for Q 2 > —m^, with the usual normal thresholds in J and u. Analyticity in t can also be established for 
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Q 2 > for the leading perturbative diagrams, and presumably holds in general. Given these results, the arguments 
of Martin fl9l - l2l| establish the Froissart bound [22j for the 7*p cross section, hence F% p - 

The phenomenological success of vector meson dominance of the electromagnetic current matrix elements between 
hadronic states provides physical evidence that the required off-shell continuation from photo-production to DIS 
electron- nucleon scattering is correct and relevant [23l425| . Furthermore, the continuations in mass- variables for 011- 
shell pion to off-shell weak-currents in PCAC-current algebra relations have had a history of remarkably successful 
predictions @, [H, HI . 

Extensive analyses of experimental data on high energy hadronic and photo-production cross sections dramatically 
demonstrate the early appearance of the In 2 s Froissart-like behavior [3, fiol lll| . This can be understood in terms 
of QCD processes at the quark-gluon level (27l - [30| . Since deep inelastic 7*p scattering is smoothly connected to ^p 
scattering by continuation in Q 2 (l8j , it is natural to assume the In 2 s Froissart behavior of the photo- production 
cross section will also appear in DIS for high 7*p energies with the substitution s — > s. In fact, detailed perturbative 
arguments [3l[ indicate that unitarity begins to be violated at remarkably small values of v = ln(l/ai), e.g. for v > 3 
at Q 2 = 10 4 GeV 2 , in the usual descri ptio n of the QCD evolution of through the Dokshitzer-Gribov-Lipitov- 
Altarclli-Parisi (DGLAP) equations |32h34| . suggesting the onset of non-perturbative Froissart-like behavior. 

The role of Mandclstam s is played in 7*p scattering by s = W 2 , the final state hadronic energy squared: 



s = W 2 = (q - 



x 



Q 2 

x) + m 2 -> — , x -> 



(1) 



For Q 2 fixed and large relative to the square of the nucleon mass, i.e. 



Q 2 > to 2 , 



Eq. (TTJ shows that a saturated In 2 



bound on the 7*p cross section translates into a In (1/x) 



bound on the small x (large 
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behavior of F^. 



In their analysis of the early ZEUS data on F% (x, Q ) [35J, [36[, Berger, Block, and Tan [1| assumed that the j*p 
cross section should show Froissart-like behavior in 1/x. The success of their model amply supports this assumption. 
In a subsequent paper Q , those authors refined their saturated "Froissart" parameterization and obtained an excellent 
global fit to both the x and Q 2 dependence of the ZEUS data, with 6 free parameters describing hundreds of points 
of data. This fit was later used along with the Feynman wee parton picture to predict the UHE v-N cross sections 
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Releasing one more parameter, the present authors lil. [il 



fit the joint ZEUS |35|,l36i and HI (37J determinations of 
the e ± p DIS cross sections as combined by those groups [17| , a combination that resolved some of the tension between 
previous individual ZEUS and HI analyses. This fit, summarized below, provides very accurate values of F2 P (x,Q ) 
over a large region of the x-Q 2 plane that includes some 335 data points. It was used in [|J in conjunction with the 
wee parton picture to predict UHE neutrino cross sections for two full quark families based on the extrapolation of 
the analytic expression for F^ v into the Q 2 > M\ and x < 10" 10 regions needed to evaluate the cross sections up to 
E v = 10 14 GeV and beyond. In our companion paper, Part II, we improve this calculation by including the b-quark 
and the NLO corrections to the wee parton limit, and extend the energy range up to the E v = I0 17 GeV energy 
reported by several neutrino detectors 114, Il5ll . 

The global fit function used in Q and jl6l|. which ensures that the saturated Froissart In 2 (1/x) behavior dominates 
at small x, takes the form 



F7M 2 



(1 



F P 



+B(Q 2 )ln 



1 — Xp 
2 fx P 1 



A(Q 2 )\n 



f Xp 1 



\ X 1 — Xp 



x 1 — xp 



where 



A(Q 2 ) 
B(Q 2 ) 



a 
b ■ 



- ax InQ 2 
&ihiQ 2 - 



-o a In 2 Q 2 , 
b 2 hi 2 Q 2 . 



(2) 



(3) 



As is evident from Eq. (p}, this form is equivalent to the quadratic expression in Ins familiar in fits to hadronic data 
with the Q 2 dependence rearranged and extended. 
At small x or large v = ln(l/a;), the expression in Eq. ([2]) becomes a quadratic polynomial in v with 



F^( Vl Q 2 ) = F^(e- V , Q 2 ) -> C Qf (Q 2 ) + C lf (Q 2 )v + C 2f (Q 2 )v 2 + 0(e~ v ). 
The coefficients Ci are again quadratics in InQ 2 , 

C* 0/ (Q 2 ) = Fp/(l-x P ) + A(Q 2 )v a + B(Q 2 )v 2 , 
C lf (Q 2 ) = A(Q 2 ) + 2B(Q 2 )v , 
C 2f (Q 2 ) = B(Q 2 ), 



(4) 

(5) 
(6) 
(7) 
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where vq = ln[xp/(l — xp)}. We will use this structure repeatedly in the analysis below. As we will see in the 
Appendix, the neglect of the terms of order e~ v in Eq. Qj, important for v ~ 0, will not affect our results at large v. 

The procedure used in fitting the combined HERA data is described in references [38| and [lH . The parameter x p 
was fixed at the value 0.11; the HERA data are sparse for larger x. Fp, the value of at xp, and the other 6 fitting 
parameters are listed in Table U together with their errors. Also shown are the renormalized minimized y 2 value (38j . 
the number of degrees of freedom and the renormalized x 2 per degree of freedom for our new analytic form for the 
combined ZEUS and HI results [n( . 



TABLE I: Results of a 7-parameter fit to the HERA combined data for F^ p (x, Q 2 ) for 0.85 < Q 2 < 3000 GeV 2 and x < 0.1. 
The Xmin is renormalized by the factor TZ = 1.1 to take into account the effects of the cut at Axf max = 6 introduced by the 
sieve algorithm used in the fit [3§ |. 



Parameters 


Values 


a 


-8.471 x 10" 


2 ±2.62 x 10 


-3 


ai 


4.190 x 10" 


2 ± 1.56 x 10 


-3 


a 2 


-3.976 x 10- 


3 ±2.13 x 10" 


-4 


bo 


1.292 x 10" 


2 ± 3.62 x 10" 


-4 


bi 


2.473 x 10" 


4 ± 2.46 x 10" 


-4 


b 2 


1.642 x 10" 


3 ±5.52 x 10" 


-5 


F P 


0.413 ± 0.003 


Xmin 


352.8 






K X Xmin 


391.4 






d.O.f. 


335 






H X Xmin/d-O.f. 


1.17 







The high quality of our Froissart-boundcd fit to data that range at the limits over ~ 5 orders of magnitude in x 
and ~ 3 orders of magnitude in Q 2 lends strong support to the proposal that the cross section for nucleon scattering 
with off-shell photons obeys the saturated Froissart bound in the 7*p Mandclstam variable s = W 2 . The errors in 
the parameters at and bi are typically a few percent except for b±; this is not well determined, with a size and error of 
the order of the errors in the other parameters. Since the fit depends linearly on the parameters, the errors propagate 
linearly, and the correlated percentage errors in the extrapolation of our fit to ultra-small x are quite small. 



III. DERIVATION OF LOW-i QUARK DISTRIBUTIONS FROM F 2 7P 
A. Relations for the quark distributions 

We start by introducing the non-singlet (NS) quark distributions [39| 

Vi = x(qi-qi), i = 1,2,..., (8) 

X3 = x{u + u — d — d), (9) 

T 8 = x(u + u + d + d-2s-2s), (10) 

T i5 = x(u + u + d + d + s + s-3c-3c), (11) 

T 24 = x(u + u + d + d + s + s + c + c- 4b -46), (12) 



and the singlet distribution 



F s = x(u + u + d + d+s + s + c + c+b + b-\ ), (13) 



where the quark distributions arc all defined at a given order in perturbativc QCD. 

We will be concerned mainly with very small x. We will take s = s, c = c, and 6 = 6, and will neglect the small 
differences between the u and d quarks at large x. The effects of the valence quark distributions u v = V\ and d v = V2 
are also quite small, and we will take d v = (l/2)u„, with u v = U, a reasonable approximation, while T3 — > (1/2)U. 
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F s is then related to the j*p structure function for different numbers nt of active quarks as 

F s (x 7 Q 2 ) = ^(x,Q 2 )-^T 8 (x,Q 2 )-^U(x,Q 2 ), n f = 3, (14) 

18 1 1 ^ 

F s (x,Q 2 ) = -F^{x,Q 2 )--T 8 {x,Q 2 ) + -T lb {x,Q 2 )--U{x,Q 2 ), nj = 4, (15) 

4^ C. C O IK 

F s (x 7 Q 2 ) = -F^(x,Q 2 )--T 8 (x,Q 2 ) + -T 15 (x 7 Q 2 )--T 24 (x,Q 2 )~-U(x,Q 2 ), n f = 5. (16) 
Here -F^cf is an expression of LO form in terms of the quark distributions, 

F^ = Y.e 2 x(q l +q l ), (17) 

i=i 

with the sum running over the active quarks. We will be concerned later with values of Q 2 above the 6-quark excitation 
threshold at m 2 , but not so large that t-quark effects are significant, so will generally take nf = 5 in the following 
discussion. 

The measured structure function is related to F%q by convolution with QCD corrections from the operator 
product expansion (40T - I421 ] . 



(x- 1 F 2 ^) + g(^e?)c 2s ® 5 , (18) 



where 1 is the unit operator and the convolution (g> is defined in Eq. (|A.2[) . The coefficient functions C 2q and C 2g are 
given in [3|| to NLO. Conversely, 



r 20 



2^ 29 



-l F gp_ ^ 
2 2tt 



The inverse operator can be evaluated using Laplace transforms as discussed in the Appendix. This requires an 
input for the gluon distribution g(x,Q 2 ). We take this from existing analyses of DIS data which give the gluon 
distribution in the HERA re gion, extending the fits to large v using quadratic extrapolations, a form implied by the 
DGLAP evolution equations (32h3^ for an input F s quadratic in v [la, HH, HH . We will assume that F^q is known. 

It is useful to note that T15 and T24 are given directly at Q 2 = m 2 , m 2 , the thresholds for c and b production where 
the c and b distributions vanish, by the F s distributions for rif = 3 and 4 at those thresholds, 

T 15 (x,m 2 c ) = F s (x, m 2 ), T 24 (x, m 2 ) = F s (x, ml). (20) 

In particular, the x dependence of the T15 and T24 is determined at the thresholds. 

We can use the expressions above to solve for the s, c, b, and light-quark distribution functions at small x where 
we can ignore valence effects and the very small splittings between the u and d distributions generated by Vi, V2, and 
T 3 and set those functions equal to zero. Then with xqi denoting the common small- a; distribution function u = d 
and with s = s, c = c, and 6 = 6, the quark distributions for n/ = 5 are 

(21) 

^Tib, (22) 
6 

^T 15 -lT 24 ; (23) 
S + ^T 16 + ir M -^). (24) 

We want to use these expressions to investigate the behavior of the quark distributions at ultra-low x and large 
Q 2 implied by our Froissart-bounded model for F™ . This requires that we know G(x,Q 2 ) = xg(x,Q 2 ) for use in 
Eq. (Ill), as well as U(x, Q 2 ), T 8 {x, Q 2 ), T 15 (x, Q 2 ), and T 24 (x, Q 2 ). We will take U(x, Q 2 ) and T s (x, Q 2 ) from existing 
parton-level fits to the HERA data, with T s (x, Q 2 ) extrapolated to small x using quadratic expressions in v = ln(l/x). 

With this input, we can determine F^ix.Q 2 ) from our fit to the HERA data using Eq. (fTT))) , and F s (x,Q 2 ) for 
rif = 3, Q 2 < m 2 from Eq. P^|) . The latter, evaluated at Q 2 = m 2 , determines Ti 5 (x,m 2 ) through Eq. ([2^1) . We 



xs 


= xqe 


"4 T - 


xc 


= xqe 


"12 T8 


xb 


= xq e 


"12 T8 
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can then calculate the evolution of T^(x, Q 2 ) to Q 2 



using the results in [45|, and then repeat the process 



to obtain T 2 ±. We will show before proceeding that the effects of non-singlet QCD evolution are small so that 
Ti5(x, Q 2 ) ~ Ti 5 (x, m 2 ) and T 2 i{x, Q 2 ) w T 24 (x, m 2 ), and will derive approximate analytic expressions for the evolved 
functions. We use the latter in our numerical calculations.. We return to the discussion of the quark distributions in 

Sec. [iv]g 



B. Non-singlet evolution 

Any of the non-singlet distributions evolves in Q 2 according to the equation (45| 

F NS {v,Q 2 )= f dwK NS (w,Q 2 )F NS (v-w,Q 2 ), 
Jo 

where v = ln(l/x), Fns(v, Q 2 ) = Fns{z~ v , Q 2 ), and Kns is the evolution kernel 

K NS (v,Q 2 ) = — / dse vs k NS (s,Q 2 ), 

27T« J-ioo+e 



kNs{s,Q 2 ) = exp 



(25) 

(26) 
(27) 



Here is the Laplace transform with respect to v of the n th order non-singlct splitting function in x, and 



T n(Q ,Qo 



4-7T 



n r Q 



Qf, 



d(ln Q u )a n s {Q") 



In leading order (n = 1), $ NS is just $/ as defined in [46[, and can be written as 



2s 



3\s+l s + 2 



y[V(.S+l)-^(l)], 



(28) 



(29) 



where V , ( z ) = r'( z )/r(z) is the digamma function. The function $^g(s) clearly vanishes at s = 0. Its only singularities 
in the complex s plane are poles at s = — 1, — 2, . . ., with the rightmost singularity at s = — 1 (where here s denotes 
a variable in Laplace space and is not the Mandelstam invariant). We can therefore move the contour of integration 
in Eq. (|2^|) to a line parallel to the imaginary axis with the real part of s to the left of s = 0, but still to the right of 
s = — 1, without encountering any singularities. 

The i ntegrand diverges for s — > — 1 and s — > 00, and has its minimum value on the real axis at the point so ~ 
near —1 where the derivative of the exponent with respect to s vanishes. Here v\ is given by 



(8/9)(7r 2 — 3)ti, so »i -> w for v large. The integrand has a value proportional to exp[wso 



exp[— v+ ^y32^YUl/3 + 0(l/vl)} at this saddle point. If we take the line of integration to run through the saddle point, 
we can estimate the integral using the method of steepest descents; the result is proportional to exp[— v+ y/32TiVi/S\. 
Given the behavior of the integrand, we see that the integral is exponentially suppressed for v 3> 32ti/3. As a 

result, the kernel K^ s is effectively zero except in a region in v extending only a distance ~ 32ri/3 from v = 0. This 
is small relative to v for the values of Q 2 and v of primary interest here. For example, T\ = 0.108 for Q 2 = 10 4 GeV 2 , 
so the constraint requires only that v 3> 1.2 or x <C 0.3. However, v > 11.5 for x < 10 -5 , the region of primary 
interest, so the integral in Eq. (|25[) samples only values of i 7 !/vs(v ~ w) very near v. 

To exploit this observation, we replace w by zero in the factor Fpf$(v — w, Qq) in Eq. (|25[) and shift the contour in 
Eq. (J2ZJI to the left of s = 0. We find that 

F NS (v,Q 2 ) w F ns (v,Q 2 ) f dwK NS { Wl Q 2 ) 

Jo 

= F NS (v,Ql) dw— dse ws k NS (s,Q 2 ) 

JO J-ioo-c 
pioo — c 



Fns(v,Q 2 )± 



(30) 
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Since Re s < 0, the first term in the integrand is exponentially small for v large and positive and can be dropped. We 
can close the contour to the right in the remaining term, and find that, for v 3> 32ti/3 and i;> 1, 

f81(v, Q 2 ) « F NS (v, Ql) e T ^>^ = F NS (v, Q%) (31) 

where the last relation uses the fact that 3?^g(0) = 0- That is, there is essentially no change in the NS distribution 
Fns{v,Qo) under LO evolution. 

This result generalizes to higher orders: $^(s) has no singularities in s to the right of s = — 1 in NLO and 
NNLO, and presumably also higher orders, and $^(0) = for all n, so, following the argument above, Fns(v, Q 2 ) ~ 
Fns(v,Qo) f° r v sufficiently large. 

The large-v part of the argument is essentially unchanged. To see that 3>^g(0) = 0, we note that $^s(s) is just 
the Laplace transform of the n th order quark splitting function, 



(*) = / dwe- w ^p-(w) = / dwe~ w ^ [P n (w) - P qq (w)j , (32) 



where we follow the notation in [39[ with the substitution of e~ w for x and P~{w) — P~{e~ w ). For s = 0, this reduces 
to 

roc 

<s(0) = J Q dwe- w [P qq (w)~P q ,(w)) 

dx (P qq {x) - P q9 (xj) = 0. (33) 

This expression vanishes as the result of quark number conservation (39j . The insensitivity of any NS distribution 
Fns(v,Q 2 ) to QCD evolution follows. We conclude that the result in Eg. (|3"Tj) continues to hold, with negligible 
corrections because of the smallness of the higher- order parameters r„ in Eq. (|2"T)l , r„ <C T\ , n > 1 . 

To see what residual effects there are in a realistic case, suppose that Fns( v , Qq) = J2 n =o c nV n > this is the asymptotic 
form of our Froissart-bounded fit to the HERA data, Eq. ^ at large v, and is also the form of any of the non-singlet 
distributions listed above. It is useful in this case to use the alternative form of Eq. (|2"5|) given by the convolution 
theorem for Laplace transforms, 

Fns(v, Q 2 ) = C- 1 [k NS (s, Q 2 )fNS,o(s); v] , (34) 

where Jns,o(s) is the Laplace transform of the initial distribution Fns(v, Qq) a t Qq with respect to v, /jvs,o(s) = 
J2n=o c n n\/s n+1 . In LO, this gives the evolved function 

FxsM*) = ^h / ^e^(«".^«. (35) 
n=0 Zm J -i°°+* s 

Since the only singularity of the integrand to the right of s = — 1 is the pole l/s n+1 at s = 0, we can shift the 
integration contour to the left of s — as shown in Fig. [IJ picking up the residue of the function exp[ws + ti^tv^ 5 )] 
at the pole, and find that 



ds n 



Ul ^ J^-e'^'&M. (36) 

-?oo — c w 



s=o 2iri 7_ < „_ / , s n+1 



The line integral which remains can be taken to run through the slightly-shifted saddle point near s = — 1 , and again 
gives a contribution which is exponentially small for v large and can be dropped. 
For an input distribution v 2 , we get an evolved distribution 



v 2 -> 



Similarly, 



r~y) Tl ] -(6+y/(l))n+0(e-*) (37) 
(u- 5.4397ti) 2 +6.8219ti. (38) 



w^w-(^-^)t 1 + OCe- 17 ) =v- 5.439771. (39) 



8 



a 




FIG. 1: Integration contours for the inverse Laplace transforms in Eq. (|35p and Eq. (|A.28[) : (a), the original contour (— ioo + 
e, ioo + e) which avoids the rightmost poles of the integrands at s = on the right; (b), the shifted contour, broken into a loop 
around the origin in the s plane, and a line integral just to the right of the singularity at s — — 1. That integral can be taken 
most efficiently to run through the saddle point near -1. There are further singularities at s — —2, . . . , indicated by dots. 



Finally, since $^5(0) = 0, a constant input function is unchanged in the evolution up to exponentially small terms. 
Combining terms and expressing the result in terms of -F/vs(i>, Qo), Eq. (|35[) gives the evolved distribution 

F NS (v,Q 2 ) = F N s{v l ,Ql)-c 2 {&+^"(l)^T 1 {Q\Ql)+0{e- v ), (40) 



8^ 2 10 , ^ ^ 2 



— -yJnW,^), (41) 

where c 2 = C2.ns{Qq) is the coefficient of v 2 in F^s( v ,Qo)- The NS evolution has simply shifted v by a small 
Ti-dcpcndcnt constant and added a small constant term. For example, v —tv'=v — 0.587 with an additive constant 
0.736 c 2 for Q 2 = 10 4 GeV 2 , Ql = 4.5 GcV 2 , n = 0.108. The constant term can be neglected for the values of v of 
primary interest here. We find, therefore, that an excellent approximation for the complete evolved distribution is 
Fns(v,Q 2 ) =F NS (v',Q 2 ). 



IV. QUARK DISTRIBUTIONS AT VERY SMALL x 
A. Results 

We will only look at the quark distributions in the region of large Q 2 , so will take n/ = 5. Then from Eq. (|16[) with 
T3 = 0, F s is given in v space by 

F. - ~ !*■ + ^ - |T 24 - ^U, n f = 5, (42) 

where Tj and U are Ti and U evaluated in v space, with x — > e~ v . 

We have determined F^q from F^ p using the NLO transformation in Eq. p^|) as described in the Appendix. The 
results are given analytically for large v in Eq. (|A.31|) . We used a gluon distribution G obtained from a fit to the CT10 
gluon distribution [43, of the form in Eq. (|4]), quadratic in v and InQ 2 . This form is necessary for consistency 
of the DGLAP evolution equations if F2(v,Q 2 ) has this form. The expression for G was fitted over the region 
2 x 10" 4 < x < 0.01 and 10 GeV 2 < Q 2 < 1000 GeV 2 , and then extended to all v, Q 2 . We note that the CT10 gluon 
distributions obtained in NLO and NNLO are very similar, and agree also with the HERAPDF results [13, S3- The 
resulting F^q, with the transformation in Eq. (fPU)) calculated in NLO, is compared with F^ in Fig. [2] The changes 
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are on the order of 5-10%, with a much smaller uncertainty from the gluon term. While we regard it as unlikely 
that higher order contributions to the functions C^q and C2 g in Eq. (fT9")l would affect the results for F^q significantly, 
we emphasize that any effects would be in the individual quark distributions, and would be insignificant for the 
relations between and the corresponding structure functions F% and FO2 for charged- and neutral-currents 
neutrino-nucleon scattering given in Eqs. (|47|) and (|51|) because of the smallness of the T' terms in those equations. 




FIG. 2: Comparison of the quark-level distribution F^q in v space (red dashed curves) with the large- v extension of our Froissart 
form fit to the HERA data, Eq. © (solid blue curves), for Q 2 = 10 4 GeV 2 (top curves) and Q 4 = 100 GeV 2 (bottom). The 
two are related by Eq. (|19[) as implemented in the Appendix through the relation in Eq. (|A.31[I . 



To get an approximate extension of the (small) function Tg to small x or large v, we use a quadratic fit to T% as 
a function of v as determined from the CT10 PDFs [U HgJ over the region 10 -5 < x < 0.003 for Q 2 = m 2 , and 
calculate it for larger Q 2 using the expression in Eq. (|40[) . We then determine T15 and T24 at the c and b thresholds 
Q 2 = to 2 ., m| using the expressions in Eq. (|20)) . again evolved to higher Q 2 using Eq. (|40[) . The relations in Eqs. (f2Tj) 
to (f2~4")l then determine the s, c, 6, and light quark distributions in terms of F^q. The results are shown in Fig. [3] for 
5 < v < 32.2, corresponding to 0.0067 >i>lx 10 -14 . For comparison, the lower limit of the HERA data is on the 
order of x = 10~ 4 , v = 9.2 for Q 2 of a few GeV 2 . 



B. Neutrino structure functions and the wee parton limit 



The assumption that the differences between the quark distributions tend toward zero for small x and large Q 2 — the 
wee parton picture — was used in 0, Q to calculate the neutrino and antineutrino cross sections on an isoscalar nucleon 
JV = (p + n)/2 at very high energies in terms of F^ as extrapolated to small x. We find here that there is not a 
proper wee parton limit. This may be seen from seen from Fig. 21 The c and b distributions remain quite different 
from the light quark distribution at small x or large v and large Q 2 , with the c/qe and b/qe ratios implied by Eqs. 
(|2"2"]) - (|24[) approaching nearly constant values at small x which decrease only logarithmically in Q 2 . 

Perhaps surprisingly, this result does not affect the supposed "wee parton" calculation of the neutrino cross sections 
in practice. The dominant structure function in charged- current neutrino scattering is F^^ is given in terms of quark 
distributions for rif — 5 by 

F 20 = x(u + d+2s + 2b + u + d+2c) (43) 

= ^ s (44) 

4^ 1 1 

= T^ F lo ™ (5^ 8 - 5T 15 + 3T 24 ) - (45) 
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FIG. 3: Plots of quark distributions determined from the Froissart-bounded fit to F^ v versus v = ln(l/a;) for Q 2 = 10 4 GeV 2 : 
top to bottom, xqt{x,Q 2 ) (green curve), xs{x,Q 2 ) (red curve), xc(x,Q 2 ) (blue curve), and xb(x,Q 2 ) (black curve). 
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FIG. 4: The ratios, top to bottom, of the distributions xs/xqe (red curve), xc/xqi (blue curve), and xb/xqi (black curve) 
plotted versus v = \n(l/x) for Q 2 = 10 4 GeV 2 . 



The results for antineutrino scattering are identical, i 7 ^ = F%. 

We obtain the observable neutrino structure function by applying the QCD corrections from the operator 



product expansion [39U42I ] to F : 



x- l F'f v) = 



1 + S C2 « 

Z7T 



-l p"W 



c 



2g<i$g, 



(46) 



where the coefficient functions C^q and C^g are given in Eq. (|A.3[) and Eq. (|A.4[) . and the effective charges C2,i are the 



11 



coefficients of the quark distributions in Eq. (|43|) . This transformation converts the F^q term in Eq. (|45l) to F^ v as in 
Eq. (|18p . The gluon term in Eq. (|4"rJ)) is absorbed in the process, and only Ci q acts on the Tj in Eq. (|4"5"1) . The result 
in v space for Q 2 > m 2 , rif = 5, is 

= (45/11)*?* - 1 (sig - 5TJ B + 3T^ 4 ) - gtK (47) 

The original functions Ti(v,Q 2 ) in u space are quadratic polynomials in u. The transformed functions Tl(v,Q 2 ) 
are again quadratics for v large. Their calculation is discussed in the Appendix, where we show that up to a small 
additive constant T'(v, Q 2 ) is simply T evaluated at a shifted value of v, 

f'(v,Q 2 ) =f{v T ,Q 2 ) + constant + O (c~ v ) , v T = v + constant. (48) 

The details are given in the Appendix around Eq. (|A.33[) . 

In the wee parton limit u = u = d = d = s = c = b = qe, the functions T[ and U' in Eq. (|47[) vanish, and F^wle = 
(45/ll)F 2 yp . The factor 45/11 = 10/ (22/9) in this expression is just the ratio of the charge factors ( J2i c 2,i) / ( Si e l) 
for rif • = 5, a relation used in the case nj = 4 in the calculations of neutrino cross section in @, As seen in Fig. [5] 

the results for F^^ obtained in this limit agree very well for large v with those calculated using Eq. (|47|) . for example, 
to ~ 3.2% (1.2%) at v = 12 (32) and Q 2 = 100 GcV 2 , with the errors decreasing with increasing Q 2 to 1.1% (0.4%) 
at v = 12 (32) for Q 2 = 10,000 GcV 2 , even though the wee limit does not really exist for the quark distributions 
derived here. These differences are just discernible in Fig. [5l and are not significant for applications at very small x 
or large v. 

It is not obvious that the relation F^^ w J^^ee = (45/11)^^ should hold as well as it does. In particular, the 
results in Fig. [4] show that the c and b PDFs are significantly smaller at all v than the light-quark PDF qg, while the 
valence distribution U = u v ps 2d v vanishes at large v. However, F^ is fixed by experiment. The overall decrease 
in the contributions of the s, c and b quarks to F^ is therefore compensated by an increase in qg. The c quark also 

appears with 4 times the weight of the s and b quarks in F^, but equal weight in F% , with the result that the 
different errors in s + 6 and c tend to cancel in the latter; somewhat accidentally, the cancellation in nearly complete. 




v = ln(l/x) 

FIG. 5: Comparison of the dominant structure function F^"^ in charged current v-N or u-N scattering calculated for, top to 
bottom, Q 2 = 10, 000, 1000, and 100 GeV 2 using the complete expression in Eq. (|470 (dashed red curves), and the approximate 
distributions F% wee w (45/ll)-F.7 P (solid blue curves) derived assuming the validity of the wee parton limit for the quark 
distributions. The limits v — 5 (32) of the range shown correspond to x = 0.007 (10~ 14 ). 

Similar results hold for the structure functions for neutral current v-N and v-N scattering. As discussed in the 
accompanying paper, Paper II, the dominant structure function is FO^ ■ The uncorrected function is given in terms 
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of the quark distributions by 

F0^ } = x(d + u + d + u) (L 2 U + Rl + L 2 + R 2 ) /4 

+x (a + b + s + b) (L 2 + R 2 d ) /2 + x (c + c) (Z 2 + R 2 U ) /2 (49) 

= Y2 {FiE ^Uj [3(L 2 + Rl) + 2(L 2 U + i? 2 )] 

~132 t 4(i ' + ^ ~ L " ~ R ^ (5Ts ~ 5TlS + 3Tm) ■ (50) 

The neutral current coefficients are L u = 1 — | sin 2 (9^, = — 1 + | sin 2 0^, i?, u = — | sin 2 and i?^ = | sin 2 
We use the value sin 2 0^=0.231 [H, [Ho| in all calculations. 

Upon transforming to the physical structure function through Eq. (|4"o]) , F^ P — > F^ v and Tj — > TV as above. The 
result in w space for Q 2 > to 2 , rij = 5 is 

F0f } = Jj (V ~ [3(LS + Rl) + 2(Ll + Rl)] 

~ W + R d) - L l - R D] ( 5f s - 5f { 5 + afu) ■ (5i) 

The naive wee parton limit, with all quark and antiquark distributions assumed to be equal, qi = q% = q, is obtained 
by taking the T/ and U' as zero, giving F0^ ee = (9/22)F 2 7P x [3{L 2 + R 2 ) + 2(L 2 + i? 2 )] . The coefficient of F 2 7P is 
just a ratio of the sums of charge factors ( c 2,i) / ( J2i e T) m Eqs. dH|) and (fTT)) . 

The exact results for FO^ m v space are compared to the wee approximation in Fig. [SJ The agreement is very 
good for v large, with agreement to ~ 3.8% (1.5%) at v = 12, (32) for Q 2 = 100 GeV 2 , decreasing with increasing 
Q 2 to 1.3% (0.5%) at v = 12 (32) for Q 2 = 10, 000 GeV 2 , even though, again, the wee limit does not really exist for 
the quark distributions derived here. These differences are quite small as seen in Fig. [51 and are not significant for 
applications to ultra high energy neutrino cross sections. 




v = ln(l/x) 

FIG. 6: Comparison of the dominant structure function FO^"' 1 in neutral current u-N or v-N scattering calculated for, top to 
bottom, Q 2 = 10, 000, 1000, and 100 GeV 2 using the complete expression in Eq. (|51|l (dashed red curves), and the approximate 
distributions FO^L ~ (9/22)F 2 7P [3(L 2 d + R d ) + 2{L 2 U + Rl)] (solid blue curves) derived assuming the validity of the wee 
parton limit for the quark distributions. The limits v = 5 (32) of the range shown correspond tox= 0.007 (10~ 14 ). 

To summarize, we have shown that the quark distributions derived from our Froissart-bounded fit to the HERA 
data do not satisfy the wee parton assumption, that all quark and antiquark distributions approach a common 
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limiting distribution at small x and large Q 2 . The quark distributions actually diverge from each other proportional 
to v 2 = ln 2 (l/x) for v large. The dominant neutrino structure functions and FO2 calculated in terms of 

the quark distributions are nevertheless given to high fractional accuracy by the "wee parton" relations F% ~ 
( Si c 2,i)/( Si e i)^2 P - This is important: F^ v is determined directly by data, and can be extrapolated to large v 
using a global fit to the v and Q 2 dependence of the data independently of the individual quark distributions. Our 
full results for UHE neutrino cross sections are discussed in detail in the accompanying paper [l3[ . 



C. Check of the wee parton limit in leading order 



In this section, we compare the results above with those obtained at ultra-small x in a conventional approach based 
on the solution of the DGLAP evolution equations [321 - I34T ] starting from input boundary functions at a hxed initial 
value of Q 2 . Our main intent is to investigate the wee parton limit from a different point of view. We first consider 
a LO calculation based on boundary functions derived directly from our Froissart bounded fit to the HERA data 
modulo small non-singlet terms taken from earlier analyses by the CTEQ and MSTW groups. Our approach follows 
[l6j |. using highly accurate and fast numerical algorithms for directly calculating inverse Laplace transforms (5ll - [53| . 
now modified to probe ultra-small x [H3 |. In particular, we can now calculate the Q 2 evolution of both singlet and 
non-singlet structure distributions down to x values as small as 10 -14 , far below the range of the HERA data, and 
determine the quark distributions needed for the calculation of neutrino cross sections at neutrino energies up to 
E v = 10 17 GeV. 

We emphasize that our intent here is only to investigate the wee parton limit. A leading-order calculation is certainly 
not reliable at the very small x and large Q 2 needed for neutrino cross sections since higher order contributions in the 
QCD coupling a s are known to be significant, and LO results also show systematic deviations from the HERA data 
[161 ] . However, if the wee parton limit fails at LO, it most likely fails at all orders. 

The essential ingredients for a solution of the singlet DGLAP equations using the procedure developed in [4f| are the 
two boundary functions F s (x, Qq) and G(x, Qq), in the range (x,l) at the initial value Q 2 = Qq. In the work in [lj| we 
constructed F s (x, Q 2 ,) from our analytic fit to the HERA data on F^, using as extra input the singlet distribution Tg 
from MSTW2008LO [H and matching smoothly to results from either CTEQ6L [H or MSTW2008LO for x > 0.03. 
Our results at small x are essentially independent of this region, and are determined by the known behavior of our 
Froissart-boundcd fit at small x. 

We chose Qq = 4.5 GeV 2 , a value in the region where the data are dense, and determined the gluon distribution 
G(x, Qq) in terms of F^ v and its derivatives with respect to InQ 2 and x by inverting the evolution equation for F^ ■ 
The details are given in [l6j . We then determined F s and the non-singlet distributions T15 and T24 by numerically 
inverting the Lap lace transforms in our analytic solution of the evolution equations in Laplace space using the tech- 
niques in [51H531 ]. and solved for the resulting quark distributions. We show the results for the MSTW-matched u, c, 
and b distributions in Fig. [7J 

The results for the quark distributions shown in Fig. [7] arc strikingly large compared to those obtained in Sec. lIV Al 
Fig. [3[ The figure, plotted on a logarithmic scale, suggests somewhat misleadingly that the naive wee parton limit, 
in which all the quark distributions approach a common limit at small x and large Q 2 with xqi(x, Q 2 ) — > xq(x, Q 2 ), 
i — u, u, . . . ,b, 6, holds in a perturbative treatment. Thus, while xc and xb differ from xu by about 4.5% and 11.5% 
at v = 18.4 (x = 10~ 8 ), the two are within a half percent and 2% of xu at v = 32.2 (x = 10~ 14 ), with the fractional 
difference still decreasing. The light quark PDFs xd, xd, xu, xu, are indistinguishable at the level of a tiny fraction of 
a percent level for small x, while xs — xs are within a percent to fraction of a percent of xu. 

However, the actual distributions do not converge to a common limit as assumed in the naive wee parton picture, 
but again diverge from each other as v 2 = ln 2 (l/x) with increasing v as required by the behavior of the initial (and 
evolved) non-singlet distributions T15 and T24. The fractional differences converge in LO only because F s grows 
asymptotically as exp [constant x \J~rv] in LO for large Q 2 and large v. The ratio Tj/F s therefore decreases sharply, 
and the differences between different quark distributions decrease strongly as fractions of the average distribution. 
This behavior is sufficient to ensure that the relation F% ~ (45/ll)F 2 yp between the electromagnetic and neutrino 
structure functions continues to be satisfied to high fractional accuracy in LO even though the wee limit is not satisfied. 

In Fig. we compare the x-dependence of F^ \x, M"z) that results from evo lving the quarks in LO from Qq — 4.5 
GeV 2 to Mf , with that found by evaluating the global fit to HERA data [l^, Eq. ©, at Q 2 = Mf over the same 
range of x. In the region of the HERA data, y < 10, Q 2 < 1000 GeV 2 , the two versions of F^ agree reasonably well, 
typically to ~ 5% as seen in Figs. 6 and 7 of [la]. However, there are systematic deviations of the LO result from the 
data which increase at large Q 2 ; a result of this trend is evident in Fig. [SI where the LO result falls noticeably below 
the fit at Q 2 = Af§, a value well above those reached in the HERA data. We note also that the % 2 of our fit to the 
data with Q 2 > 2.7 GeV 2 is substantially better than that of the best LO result, 295 for 296 degrees of freedom as 



14 




18 20 22 24 26 28 30 

v=ln(l/x) 



FIG. 7: The ultra low x behavior of the LO quark distributions xu, xc and xb at Q 2 = M\ /3 for boundary functions derived 
from the Froissart-bounded fit to F^ v at Qq = 4.5 GeV 2 as described in [1(1 ■ Top to bottom at v = 18.4 (x = 10 -8 ), xu, xc 
and x6. 
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FIG. 8: The comparison of the results for p2 P (x,M%) obtained using the LO BDHM quarks of [13] (solid blue curve), and 
those calculated directly from the Froissart-bounded fit function, Eq. ([2]) (dashed red curve), as a function of v — ln(l/ai). The 
interval shown, 3 < v < 32.2 corresponds to 0.05 > x > 10 -14 . 



compared to 502 (1480) for the LO fit matched to MSTW2008 (CTEQ6L) at large x, or x 2 /d-o.f. « 1 compared to 
1.7(5). 

The evolved LO distribution for is clearly not consistent with the expected Froissart-bounded behavior at large 
v. While NLO corrections are known to be large and may damp the growth of the gluon distribution and Fj 9 at 
small x, the large deviation of the LO results from those given by our Froissart-bounded fit to the data suggests that a 
perturbative description may well fail, and that nonperturbative effects are likely to become important at ultra-small 
x. 
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D. Uncertainties in extrapolations of quark distributions 

A different source of uncertainty in standard analyses which affects predictions for neutrino cross sections at very 
high energies is the use of arbitrary, often power-law based, parametrizations of input PDFs that cannot be extended 
reliably to small x, a problem we have avoided by working directly with F^ p , which we assume can be extrapolated 
reasonably to small x using a theoretically consistent, Froissart-bounded form. To illustrate the problems that can 
be caused by the extension of typical parametrizations outside the region of the HERA data, we compare the present 
results with those of a LO pQCD analysis based on the parametrizations used in the comprehensive MSTW paper 
[55| . We use their MSTW2008LO parton parameterizations with the published parameter values to determine their 
boundary functions at Qq = 1 GeV 2 . MSTW numerically solved the coupled DGLAP equations on a fixed grid, 
1CP 6 < x < 1, 1 < Q 2 < 10, 000 GeV 2 , and fixed the parameters by matching the results to data. 

We are primarily interested in the behavior of PDFs parametrized this way well below their grid limit on x, x > 10~ 6 . 
We assume that the same ad hoc starting expressions, satisfactory on their grid, can be extended to x = 10~ 14 region, 
an obvious source of potential problems. For example, the MSTW2008LO fit function for the gluon PDF uses a single 
inverse power of x, xg(x,QQ) oc a: -0 ' 837 , to handle the small x description, a Froissart-violating form that dominates 
the evolution at small x and drives the behavior of the quark distributions (57j . Such assumed power-law dominated 
forms for the gluon distribution are common in the literature. 

We solve directly for the quark PDFs using the same methods as above, but in terms of the MSTW2008LO input 
functions instead of those derived from our fit to '. In Fig. [91 we display the resulting small x behavior of the 
up, charm and bottom quark distributions at Q 2 = M| over the range from x = 10~ 14 up to x = 10~ 8 . The 
u,u,d,d, s and s quarks are indistinguishable on the scale of the figure. For fixed Q 2 , the quark PDFs rise as an 
inverse power law xq oc a; -0 ' 81 as x decreases, are many orders of magnitude larger at small x than those obtained in 
our earlier analysis, Fig.[3J and in our LO analysis starting with boundary functions determined by F^ 9 at Qg = 4.5 
GeV 2 , Fig. [7] They are presumably unrealistically large. 

A second remarkable feature of Fig. [5] is that the c and b quarks never converge to the others, all rising with the 
same approximate power law, but with different normalizations at large v. The quark distributions are not consistent 
with the wee parton picture, since the ratios of different quark distributions arc fixed, and the different distributions 
do not converge toward a common value but actually diverge from each other as a power of l/x as x — > 0. 
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v=ln(l/x) 

FIG. 9: Leading-order solutions for ux, xc, and xb at Q 2 = M§/3 as functions of v — ln(l/ir) for initial gluon and quark 
distributions defined by extending the MSTW2008LO parametrizations, as given, from the HERA region to the ultra-small x 
region below x = 10 -8 , v > 18.4. Top to bottom, ux, xc, and xb. 

The power-law growth of the quark distributions at small x is what is expected for an a; -0 84 power-law input for 
the gluon since the gluon drives the small x behavior of the quark distributions [HtJ • The divergence in G at small 
x is sufficiently strong that it already dominates the very small x behavior of F s at Q 2 — m 2 , the threshold for 
charm production where F s determines the initial non-singlet distribution T 15 , T 15 (x, m 2 ) = F s (x,m 2 ). As a result, 
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the gluon power divergence at small x is passed on to the difference between the charm and light quark PDFs even 
though there is no direct coupling of T 15 to the gluon. The same is true, with increased effect, at the b threshold, so 
xb(x, Q 2 ) also rises with approximately the same power for x —> as is evident in Fig. H3 and the ratios of the quark 
PDFs do not converge as required in the wee limit. 

As a result of this universal behavior, F^ 9 and the neutrino-nucleon structure function show the same gluon- 

driven power law growth at small x. In terms of our discussion in Sec. QJ the virtual gauge boson-nucleon cross 
sections rise with s as s a with a ~ 0.8, in radical violation of the Froissart bound [l9l - l22| . a behavior which illustrates 
the danger in extending ad hoc paramctrizations outside the region in which they are constrained by data. 

A related problem appears in the extrapolation of quark PDFs from the experimental region to very small x when 
this is done by matching assumed forms without a firm theoretical basis to the fitted PDFs near the boundary of 
the experimental region. For example, the form a + b\n(l/x) used to extrapolate the published NLO MSTW quark 
distributions from the limiting grid value of 10~ 6 in MSTW2008 to the values ~ 10~ 10 needed in the calculation of 
neutrino cross sections in |58j], is inconsistent with the power law rise of the actual LO and NLO "MSTW" quark PDFs 
at ultra-small x. In the LO case checked here, the function used leads to a major underestimate of the extrapolated LO 
distributions and the corresponding high energy cross sections relative to those obtained in the direct (but unrealistic) 
LO calculation. Interestingly, the extrapolation of the NLO PDFs in [f38| still gives a cross section larger than that 
obtained using our Froissart-bounded extrapolation of F% ', quadratic rather than linear in ln(l/cc), and the relation 

F% « (45/ll)F^ rp ; this is presumably the result of matching to the rapid growth of the power-law variation of the 
MSTW quarks, and the limited range of the extrapolation. 

V. CONCLUSIONS 

We have investigated the results that follow from the assumption that the cross section for the scattering of a 
virtual photon with q 2 = —Q 2 < from a nucleon, hence the structure function F^ v in deep inelastic e p scattering, 
is hadronic in nature with the same Froissart-bounded structure as is observed in hadronic and real ^yp scattering. 
We have presented theoretical arguments in favor of this assumption, which is supported experimentally by the very 
accurate fit to the HERA data on F^ 9 obtained in earlier work. This fit is quadratic in the natural variable v = ln(l/x) 
for x small, and allows a reliable extrapolation of F^ 9 to ultra low values of x. 

We have used this fit in conjunction with information on the small non-singlet function Tg and the gluon distribution 
extrapolated consistently from results of the CT10 analysis of the HERA data [ijj to derive a complete set of quark 
distributions for n/ = 5 active quarks for x > 10 -14 or v < 32. The derivation does not use the DGLAP equations, 
which are expected to break down at very small x [3l| . These quark distributions do not show the limiting behavior 
expected in the wee parton picture, in which the deviations of the distributions from one another tend to zero at small 
x and large Q 2 , but actually diverge for x —> 0. 

We show that, despite the failure of the wee parton picture at the quark level, the relations between F^ 9 and the 

dominant structure functions F^ and FQ 1 ^ in charged- and neutral-current neutrino scattering derived in the wee 
parton picture continue to hold to high accuracy at very small x. With this established, the use of the (supposed) 

wee parton relations to predict the dominant neutrino structure function F%^ in terms of F^ 9 gives results that 
are independent of assumptions about the gluon distribution, and appear to hold quite generally. These relations are 
exploited in the accompanying paper [l3|. The neutrino cross sections will be accessible at energies E v up to 10 17 GeV 
in planned neutrino observatories, requiring values of x down to x = 10~ 14 . This corresponds to a relatively modest 
extrapolation by a factor of ~ 3 in v from the upper values v ~ 10 explored at HERA to the values v ~ 30 needed 
for E v ~ 10 17 GeV. We emphasize that, through the connections established here, measurements of the neutrino 
cross sections would allow, through the structure functions, the exploration of hadronic interactions at energies not 
otherwise accessible. 

As part of our analysis, we obtain simple analytic expressions for the effects of non-singlet DGLAP evolution on 
the functions Tg, T15, and T24 needed in the derivation of quark distributions, and of the effects of the NLO QCD 
corrections needed to transform between the LO-type expressions for F^q in terms of quark distributions, and the 
physical F^ 9 . 

We also compared our general results to those obtained in a LO analysis based on boundary functions obtained 
up to the small non-singlet term Tg entirely from our Froissart-bounded fit to the HERA data. The LO results are 
markedly larger at very small x and large Q 2 , and are not consistent with the extrapolated fit. There is again no 
proper wee parton limit for the quark distributions. 

Finally, to illustrate the difficulties that can ensue from the extrapolation of the ad hoc parametrizations of parton 
distributions used in standard analyses of deep inelastic scattering out of the experimental domain where the param- 
eters are determined, we calculated the "MSTW" quark distributions which follow from a direct extension to very 
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small x of the analytic expressions used to specify the initial boundary distributions in the MSTW2008LO fit to F^ v ■ 
The results are drastically different from the Froissart-bounded results found here, and also from the a + b\n(l/x) 
extrapolations of the fitted distributions MSTW distributions used in the prediction of neutrino cross sections in [58| . 
It appears to be far better to extrapolate directly, and use the results to determine boundary functions for a 
DGLAP analysis, or, as here, to determine the PDFs from the extrapolated global fit to the v- and <3 2 -dependcnt 
data. 

We conclude that the cross sections and quark distributions calculated from the small x, large Q 2 extrapolation 
of F^ p {x : Q 2 ) from our saturated Froissart-bounded fit to the HERA data are the most physically motivated, are 
consistent with all other hadronic cross sections including jp, and provide the best estimate of the UHE energy 
neutrino-nucleon cross sections, which we develop fully in Part II [l3j . 



Acknowledgments 

M. M. B. and L. D. would like to thank the Aspen Center for Physics, where this work was supported in part by NSF 
Grant No. 1066293, for its hospitality. M. M. B. would like to thank Prof. A. Vainshtein for valuable discussions. P. 
H. would like to thank Towson University Fisher College of Science and Mathematics for support. D. W. M. receives 
support from DOE Grant No. DE-FG02-04ER41308. 



Appendix: Calculation of F^'q 



We recall that the structure function is given in terms of the quark-level expression F^q in Eq. (|T7|) by convo- 
lution with a set of coefficient functions from the operator product expansion [4ll [42j , 



x-^r 



1 + ^ 

Z7T 



(*- 1 ^) + £(E«?)<v 



9, 



(A.l) 



where the convolution (g> of operators A and B is defined as 

- 1 dz 



A®B = 



-A(x/z)B(z) = 



1 dz 



A(z)B(x/z). 



(A.2) 



The operator 1 in Eq. (|A. 1|) is the unit operator and the sum over charges in the second term runs over active quarks 
and antiquarks. 

The coefficient functions Ci q and Ci g depend on the renomalization scheme used in perturbative calculations and 
the order to which they are carried. We assume the use of the standard MS scheme in which, at NLO [39| , 
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/ln(l - z) 



1 - z 



C 



+3 + 2z- (l + z)ln(l-z) 

1-2 



1-z 



■ In z 



2;/ 



((l-z) 2 + z 2 )ln- 



8: 



8z - 1 



(A.3) 
(A.4) 



The coupling a s (Q 2 ) is to be evaluated at the same order. 

The expression in Eq. (|A.1[) is usually used to determine a; _1 -F^ p from the individual quark and gluon distributions 
found in fits to the DIS data. However, the relation can also be inverted to determine x^F^q directly at a given 
order in a s in terms of the observable structure function .t _1 F s 7 p and a given gluon distribution g(x, Q 2 ), i.e., 



x r 20 ~ 



<E> x~ L F. 



G 



(A.5) 



This is the result we need to obtain the singlet quark distribution F s and individual quark distributions as outlined 
in Sec. MI Al As discussed there, F s is determined (except at very low Q 2 ) by F^q and the non-singlet functions T15 
and T24, themselves related to F s . 
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Wc sketch here the evaluation of the inverse operator and the final expression in Eq. (|A.5[) using Laplace transform. 
This requires several steps. We first multiply by x and use the second form of Eq. (|A.2[) to recast Eq. (|A.1[) in the 
form 



F^(x,Q 2 ) = [l + ^(zC 2fl )l®i&> + g + (X;e?)(*a 



2- 



G 



C 2q (z)F%(x/z,Q 2 ) + (j2ef)zC 2g (z)G(x/z,Q 2 ) 



(A.6) 
(A.7) 



where G(x,Q 2 ) = xg(x,Q 2 ). Wc next transform the terms in Eq. (IA.3I) that involve distributions. Since F{x/z) 
F^ix/z, Q 2 ) = for z < x, we may take the lower limit of integration in the convolution as 0, and write 



dzF{x/z) 
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\z F ^ 



(1-*)h 



1 F(x/z) - F(x) 
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F[x) In h x 



1 - z 

1 (F{y) F(x)\ dy 
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dw In (l - e- (v - w) ^ 



y x J y-x 

dF(w) 
dw ' 



(A.8) 
(A.9) 
(A.10) 
(A.ll) 
(A.12) 



where we have used the definition of the "+" operation in Eq. ()A.9|) . evaluated the integral on the interval (0,x), 
changed the integration variable z to y — x/z in Eq. (jA.llj) . and finally introduced the natural variables v = ln(l/x) 
and w = ln(l/y) and integrated by parts in Eq. (|A?12|) using a limiting procedure as sketched in [4(|. The function 
F(w) is defined as F(w) = F$(w, Q 2 ) = F^ p (e~ w ,Q 2 ). 

A similar calculation for the term in Eq. (|A.3|) proportional to (ln(l — z)/(l — z)), gives 



dz F{x/ z) 



/ln(l - z) 



dwln 2 { l- e -("-™) 



dF(w) 
dw 



(A.13) 



+ 

Using these results and transforming the remaining terms in Eq. (|A.7|) to v space, we obtain the expression 



F 2 ->,Q 2 ) = F%{v,Q 2 ) + ^ 

Z7T 
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6 + 7 y-) F^{v, Q 2 ) + dwH q (v~ w)F%(w, Q 2 ) 



dw 

as(Q 2 ) 
2tt 



- In 2 ( 1 - e - {v - w) ) - 4 In [ 1 - e -( v - w *> 



dF^(w,Q 2 ) 
dw 



dwH g (v-w)G(w,Q 2 ). 



(A.14) 



Here G(w, Q ) = G(e W ,Q 2 ). The functions H q and H g are defined as 

H q (v) = e- v C' 2q (e- v ), (A.15) 
H g (v) = e- v C 2g {e~ v ), (A.16) 
where C 2q (z) contains the terms in C 2q other than the delta function and the "+" terms treated above, i.e., 



CM = - 



1 + z 2 

'i + 2z In z - (1 + z) ln(l - z) 

1 — z 



(A.17) 



The right-hand side of Eq. (|A.14j) is a sum of convolutions in v space, and can be factored by Laplace transformation 
into a sum of the products of the transforms of the functions in those convolutions, 



his) - Ms) + g/ 20 (s) f-6- \f + h gl (s) + s h q2 (s)) +^±(Y / e 2 )g(s)h g ( s 



2tt 



(A.18) 
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Here f 2 o(s), f2(s), and g(s) are the Laplace transforms of F2O) F2, and G with respect to v, with their Q 2 dependence 
suppressed, 



while 



h q i(s) = C 



sh q2 (s) 



h g (s) = 



4 
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sC 



/ 20 (S) 



H q (v);s 
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(A.19) 
(A.20) 
(A.21) 
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+ C(2,s + l) + C(2,s + 3) 
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H g {v);s 
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4 ff s+2 + 4 
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(A.22) 



(A.23) 



(A.24) 



In these expressions, _ff s = ^(s + 1) — ^(1), ^>(s) = F'(s)/r(s), and £(2, s) = X^o(^ + s ) 2 ^ s * ne Hurwitz generalized 
zeta function of degree 2. The factor s which multiplies h q 2(s) in Eq. (| A. 18|) and Eq. (|A.23[) arises from the derivative 
of F 2 7 (f in Eq. (|A~T4| and th e relat ion C[df{w)/dw; s] = s C[f{w); s}. 

In the expression in Eq. (|A.18[) . f-2(s) is known from our fit to the HERA data, and g(s) is assumed also to be 
known, for example, from the extension of G(x, Q 2 ) from earlier parton level fits to the data as extended to small x. 
Solving for f2o(s), we find that 



/ 20 (S) 



/2(s)-^(E e 0^ S ^ 



[l + (o,/27r)d(«)], 
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d(s) = -6 - -7r 2 + h ql (s) + sh q2 (s). 
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Thus, inverting the Laplace transform in Eq. (|A.19|) . we find that 

■M 

f2(s) 



(A.25) 
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1 + (a s /2ir)d(s) 



2tt 
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h g (s)g(s) 



1 + (a s /2ir)d(s) ' 



(A.27) 



The inverse Laplace transform in Eq. (|A.27|) can be calculated simply analytically for v large or x small. In 
particular, in our Froissart bounded model, '(v,Q 2 ) and G(v,Q 2 ) are essentially quadratic polynomials in v for 
v >> 1 as in Eq. @. In the polynomial terms, v n — > n\/s n+1 under Laplace transformation. The exponentially small 
terms omitted in Eq. @ give extra poles for s — > —1, —2, . . ., and, if retained, lead only to terms of order e~ v or 
smaller in the final result. The main contributions for rif = 5 and v >> 1 are therefore given by integrals of the form 



2ttI 
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— 100-j-e. 



ds_ eVS [l, (a s /27r)(22/9)h g (s)} 



„n+l 



1 + {a s /2n)d(s) 



n = 0,1,2, 



(A.28) 



where the numerator function in Eq. (|A.28|) is 1 for the f 2 term in Eq. (|A.27|) and (a s /27r)(22/9)/i g (s) for the g term. 

The numerator functions have no singularities in the complex plane to the right of s = —1. The function d(s) 
has second-order poles for s — > —1, —2, —3, . . ., but these cause no problems. However, the complete denominator 
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function has a pair of complex conjugate zeros near and to the right of —1 where 
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(A.29) 



(A.30) 



e.g., at s = -0.9962 ± 0.1750 i for a s (M|) = 0.118. There are further pairs of conjugate poles near s = — 2, — 3, 

These pole positions will be shifted slightly and new poles introduced when the coefficient functions C 2q and C 2g are 
evaluated to higher orders in a s , introducing higher order contributions in l/(s + 1), but the rightmost singularities 
from the (generalized) factor 1/ [1 + (a s /2n)d(s)\ should remain very close to s = — 1. 

We conclude that the contours of integration in Eq. (|A.28|) can be shifted to the left in the complex s plane as in 
Fig.[TJto run through saddle points close to s = —1, but just to the right of the complex conjugate poles in Eq. (|A.30|i . 
picking up the residues of the integrands at s — and leaving a residual integral which is suppressed by a factor 
k e"", very small for v large. We drop the latter. 

The calculation of the residues of the poles at s = is straightforward. Thus, for the "1" term in Eq. (|A.28[) . the 
quadratic form of our input function i<2 (i>) in v is reproduced in F^qIv) with small shifts in v and an added constant 
in the v 2 term as given below in Eqs. (|A.31[) and (|A.32|) . The analytic forms of the coefficients in these expressions 
are known, but are too complicated to record here. Similar results hold for the "g" term. 

Combining the results, we find that, for general values of b(Q 2 ) — a s {Q 2 )/2-K 1 



FHo(v, Q) = F™{v f ,Q 2 )- (5.523- 17.665 6) bC 2f {Q 2 ) 

-^bG(v g ,Q 2 ) + ^b (7.549 + 5.523 b - 17.665 b 2 ) C 2g (Q 2 ) + O (e 



(A.31) 



to NLO, where C 2 f(Q 2 ) is the coefficient of v 2 in F^ p (v,Q 2 ), Eqs. (j4} and (0, and C 2g (Q 2 ) is the coefficient of the 
corresponding term in G(v,Q 2 ). The shifted arguments Vf and v g are 



<7 



4.203 6, 



= v - 4.623- 4.203 6. 



(A.32) 



The main uncertainty in the overall result for F 2 q arises from the uncertainty in the gluon distribution. This was 
treated as in Sec. IIVI using an extrapolation of the CT10 [13] G(x,Q 2 ) quadratic in v, with coefficients quadratic in 
InQ 2 , fitted to the NNLO G over the region 2 x 10" 4 < x < 0.01, 10 GeV 2 < Q 2 < 1000 GeV 2 . As noted earlier, this 
agrees very well with the HERAPDF version of G. 

We emphasize that the form of these results, with quadratics in v transformed to quadratics up to exponentially 
small corrections, is quite general, the result simply of the calculation of residues at s = 0, with all other singularities 
of the integrands, from either the kernel functions in Laplace space or the forms of F^ or G for v ~ 0, displaced at 
least to the vicinity of s = — 1. 

The calculation of the neutrino structure functions F^^ and FQ^' also requires the evaluation of the action 
of [1 + (a s /2ir)C 2q \ on the functions Tg, T15, and T24. In v space, these are quadratics in v for v large. The 
resulting transformation of the powers v n , n = 0, 1, 2 is just the inverse of that associated with the "1" term in the 
transformation F^ —J- F^q discussed above; G does not enter. Thus, 



fi(v, Q 2 ) -> f z (v T , Q 2 ) + (5.523 - 17.665 6) 6C , 2jTl (Q 2 ), 



Vt 



4.203 6, 



(A.33) 



with C 2} Ti(Q 2 ) the coefficient of the quadratic term in v in Ti(v, Q 2 )- The functions Ti{v, Q 2 ) are given in terms of 



the initial distributions Tj(u, Qq) determined at Qq 



by the expression in Eqs. (J25J), (gDJ), and (|4"Ijl . 



The calculation of the complete neutrino cross sections also requires the structure functions xF 3 and Fl ■ These are 



given to NLO, using the form analogous to that for F^ p in Eq. (|A.6|) . by 
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where, for example, Fg = u + d + 2s + 2b — u — d— 2c for 71/ = 5. The coefficient functions are 

C 3q (z) = C 2q (z)-^(l + z), C 3g = (AM) 

C Lq (z) = |z, C Lg {z) = 2z{l-z). (A.37) 

Transforming Eq. (|A.34[) to v space and factoring the resulting convolution with a Laplace transform, we find that 

h = [1 + bd(s) + bh 3q (s)] ho- (A.38) 
where / 3 and /30 are the Laplace transforms of F 3 and F 30 with respect to v, and 

*■<■>— 5 (jtt+jts)- (A ' 39 » 

Since F 3 q is a quadratic in v for v large, we can calculate the inverse Laplace transform of f 3 as above by calculating 
the residues of the integrand e vs [1 + bd(s) + bh 3q (s)] (n!/s™ +1 ) for n = 0, 1, 2, corresponding to inputs v n . The 
results give 

F 3 {v,Q 2 ) = (1 - 2b)F 30 {v 3l g 2 ) + (2 - 523 ~_ 3 ^ 499b)b C 2 , 3 (Q 2 ) + Q (e~") , (A.40) 
5 - 870& 

W3 = W+ T^26' ( ^ 

to NLO, with (\ 3 (Q 2 ) the coefficient of w 2 in F 30 (w, Q 2 ). 
Similarly, for i^, we find that 

h = bh Lq (s)f 20 + 2n f bh Lg (s)g, (A.42) 

hi « = 57T^' ^ = 7T2-7T3' {A - 43) 

where 5 is the Laplace transform of G(i>, Q 2 ). Using the quadratic forms of /20 and g in v and evaluating the residues 
at s — in the inverse Laplace transform, we get 

in NLO, with C , 2 , /0 (Q 2 ) and C 2g (Q 2 ) the coefficients of w 2 in F 20 {v,Q 2 ) and G(w,Q 2 ). 
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